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ABSTRACT 


The  fundamental  equations  of  the  transonic  box 
method  were  derived,  based  on  the  representation  of 
the  velocity  potential  by  a  doublet  distribution.  They 
form  the  basis  of  a  systematic  method  of  treating  an 
oscillating  wing  at  M  =  1,  analogous  to  the  supersonic 
Mach  box  method. 

A  digital  computer  program,  written  in  Fortran  IV, 
is  presented.  The  program  applies  to  a  planar  wing  of 
polygonal  planform,  with  a  straight  trailing  edge,  and 
as  many  as  three  st/eep  angles  along  the  leading  edge. 
For  a  maximum  of  ten  modes  of  oscillation,  the  pro- 
gram  computes  the  oscillatory  potentials  and  pressures 
and  a  generalized  force  matrix. 

Results  obtained  from  the  program  are  compared 
with  existing  theoretical  and  experimental  values. 
Several  possible  extensions  of  the  method  are  described. 
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BY 
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d 

DA 

f 
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Local  speed  of  sound;  speed  of  sound  at  infinity 
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Part  of  AY  (J) 

Pressure  coefficient 

Column  matrices  used  in  least  squares  surface  fits 
Dimensionless  length  of  box  side 
Coefficient  in  deflection 


The  data  array 

Function  which  describes  the  wing  deflection 
Factor  which  gives  y  the  proper  edge  behavior 
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surfaces 
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k 
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Indexes  specifying  box  position 
Indexes 
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Ug,  Air  speed  of  the  wing;  speed  of  flow  at  infinity 
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1.  INTRODUCTION 


The  transonic  box  program  is  designed  to  calculate  the  unsteady  potentials 
for  a  given  set  of  modep  of  wing  oscillation  and  to  compute  the  generalised 
forces.  Pressure  distributions  may  be  obtained  from  the  potentials. 

A  planar  wing  with  a  straight  trailing  edge  is  assumed.  The  oscilla¬ 
tions  are  assumed  to  be  symmetric  in  the  spanwise  coordinate  y.  None  of 
these  assumptions  is  necessary  for  the  method.  (See  Section  5.  ) 

The  basic  step  in  the  box  method  is  the  solution  oi  the  system  of 
simultaneous  equations  [Equation  (33)  ]  which  determine  a  set  of  values  of 
potential  on  the  wing  from  a  corresponding  array  of  upwash  values.  A  sur¬ 
face  is  fitted  to  these  values,  giving  a  functional  representation  of  the 
potential  that  is  used  subsequently  to  find  pressures  and  generalized  forces. 

The  method  used  is  suggested  by  the  success  of  supersonic  box 
methods  (References  1  through  4).  The  potential  is  generated  by  a  do  iblet 
distribution  rather  than  by  a  source  distribution  because  the  latter  method 
would  involve  diaphragm  regions  of  infinite  extent,  whereas  the  doublet 
distribution  is  confined  to  the  wing  and  its  wake. 


Muuncfipt  re  lea  ted  by  the  author*  15  September  1964  for  publication  at  an  RTO  Technical  Documentary 
Report. 
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2. 


THEORETIC. ..  OEVELOPMENT  OF  THE  METHOD 


1.  THE  DIFFERENTIAL  EQUATION 


We  consider  an  oscillating  body  moving  at  speed  U®  through  a 
nonviscous  fluid.  From  the  point  of  view  of  a  moving  coordinate  system 
(x,y,  z)  in  which  the  average  position  of  the  body  is  fixed,  there  is  a  flow 
past  the  body  with  velocity  U®  at  infinity.  Assume  that  the  flow  is  irrota- 
tional;  then  the  velocity  field  of  the  flow  is  the  gradient  of  a  potential  function 
$,  which  satisfies  the  differential  equation 


V2  « 


*tt  +  2  ‘  V*t  + 


V  )  1/2  (V*) 


=  0 


(1) 


(See  Reference  5,  p.  193'  where  a  is  the  local  speed  of  sound. 

Suppose  that  the  flow*  is  approximately  uniform  in  the  direction  of  the 
positive  x-axis.  This  may  be  true,  for  example,  if  the  body  is  almost 
plane  and  the  oscillations  are  small.  Then  #  may  be  broken  up  into  several 
parts.-  as 


*  =  Ucfx  ♦  f  +  ?)  (2) 

where  the  first  term  gives  a  uniform  flow,  the  second  term  gives  the  cor¬ 
rection  for  a  steady  flow  about  the  body,  the  third  term  gives  the  correction 
to  this  for  the  oscillating  body,  and  f  and  ?  are  smalL 

To  the  first  order,  f  and  f  are  different  solutions  of  the  same  differ¬ 
ential  equation 


(3) 


where  M.a  are  the  Mach  number  and  speed  of  sound  at  infinity.  (See  Refer¬ 
ence  5.  p.  198.  )  f  is  a  periodic  function  of  t.  Since  the  differential  equation 
is  linear,  we  may  put  +  -  (x,  y,  s)  elw*,  where  w  is  the  angular  frequency 

of  oscillation.  In  terms  of  the  nondimensional  quantities. 
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X 


=  x/b 
y  =  y/b 

z  =  z/b 

k  =  db/U 

CD 

(b  is  a  characteristic  length  of  the  body);  Equation  (3)  becomes 

a-M2)*  +-5  -  2iM2kp  +  M2kS  =  0  (4) 

xx  yy  zz  x 

For  M  =  1,  this  reduces  to 

p  +  p  -  2  iky  +  k2p  =  0  (5) 

yy  zz  x 

the  linearized  transonic  equation  (see  Reference  6,  p.  7 X  It  has  been  sug¬ 
gested  by  Landahl  (Reference  6}  that  the  proper  equation  to  use  instead  of  (4) 
is 

p  +  p  -  2  i  M^kp  +  M2k2  p  =  0 

yy  zz  x 

if  k»  |M-  1 1.  Comparison  of  this  equation  with  (5)  leads  to  a  similarity  rule 
for  flows  in  the  transonic  range  (see  Reference  6,  p.  18). 

The  range  of  validity  of  this  equation  is  duscussed  by  Landahl  (Refer¬ 
ence  6,  Chapter  1).  First,  there  is  the  requirement  for  linearization  in  any 
speed  range,  that  the  perturbation  potential  *  be  small.  This  is  not 
satisfied  at  the  leading  edge  of  a  wing  for  any  realistic  cross-sectional 
shape;  however,  it  may  be  satisfied  over  the  rest  of  the  wing,  if  the  wing 
has  small  thickness,  and  the  results  on  parts  of  the  wing  away  from  the  lead¬ 
ing  edge  are  not  much  affected  by  the  error  there. 

Another  restriction  peculiar  to  transonic  speeds  is  associated  with  the 
absence  of  the  term  inp^.  The  actual  flow  has  some  variation  in  local  Mach 
number  which  may  influence  the  nature  of  the  flow  considerably  if  M  is  near  1. 
The  presence  of  the  term  in«Fx  tends  to  reduce  this  influence,  but  for  k  small 
or  zero,  the  equation  is  valid  only  for  a  highly  swept  wing  with  a  pointed  nose. 

The  difference  of  the  local  Mach  number  from  the  value  1  assumed  in 
Equation  (5)  may  come  from  two  sources:  (1)  wing  thickness,  and  (2)  a  change 
in  the  free  stream  Mach  number.  Thus,  for  any  value  of  k,  there  are  limits 
on  the  thickness  ratio  and  the  Mach  number  range,  which  increase  with  k. 
Estimates  of  these  limits  are  not  possible,  because  of  the  small  amount  of 
experimental  data  available. 
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2.  BOUNDARY  CONDITION* 


v.  --- *r-  ~  %? 


The  solution  of  Equation  (1)  must  give  a  velocity  field  which  is  such 
that  a  particle  at  the  body  surface  moves  along  the  moving  surface.  If  the 
equation  of  the  surface  is 

S(x,y,*,t)  =  0 

this  equation  must  be  satisfied  when  (x,  y.  z)  moves  with  the  velocity 
Differentiating  with  respect  to  t  gives  the  condition 

ftC 

V®  -  V  S  +  ~  =  0  (6) 

ot 

This  determines  the  normal  velocity  at  the  surface. 

Now  suppose  that  the  body  (to  be  referred  to  henceforth  as  a  wing)  is 
almost  planar,  lying  almost  in  the  xy- plane  (see  Figure  IX  For  vertical 


Figure  1.  Coordinate  Systems 

oscillations  of  the  body,  rhe  upper  and  lower  surfaces  may  be  represented 
by  the  equations 


=  gu  <*.y> + 


iut  -  -  - 
e  f  (x.y) 
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■iiilUiiia kaitev»;>miU^uiak4 .  •-  >|  <**»«  *  .it r.^y.-  aamwi^Byi i i^wri ri t'ttSF* » 


Z 


=  (x,  y)  +  elwt  f  (x,  y) 

where  the  functions  gu.g£  are  associated  with  the  deviation  of  the  shape  of 
the  body  from  planar,  and  f  depends  on  the  mode  of  oscillation.  Then  on 
the  two  surfaces,  we  may  take 


Use  these  expressions  for  S  and  Equation  (2)  in  Equation  (6).  Neglect 
ing  terms  that  involve  products  of  *  or  $  with  gn,gf,  or  f,  the  resulting 
equation  may  be  broken  up  into  a  steady  part,  which  gives  the  boundary 
condition  for  f,  and  an  unsteady  part,  which  gives  the  boundary  condition 
for  The  unsteady  part  is 


ay  _  af 

d  z  dx 


+  ikf 


(7) 


where  f  =  f/b.  To  the  present  degree  of  approximation,  this  condition  should 
be  applied  at  z  =  0,  over  the  region  of  the  xy- plane  on  which  the  body 
projects. 

3.  THE  BOUNDARY  VALUE  PROBLEM  FOR  * 


In  linearized  theory,  a  disturbance  of  a  flow  at  Mach  1  does  not  have 
any  influence  upstream.  Consequently, 

■?(x,  y,  z)  =  0,  x  <  0  (8) 


if  the  body  lies  in  the  region  x  2  0.  This  is  one  of  the  conditions  y  must 
satisfy. 


?  is  a  solution  of  Equation  (5)  in  all  space  outside  S  and  W,  the  regions 
in  the  xy-plane  occupied  by  the  wing  and  its  wake  (see  Figure  1).  In  general,  j 
y  is  discontinuous  in  these  regions.  A  boundary  condition  on  W  is  obtained 
by  equating  the  pressures  above  and  below  the  surface  of  the  wake.  From 
the  linearized  form  of  the  pressure  coefficient, 

l 

C  =  "2  »  +  iky)  f 

P  *  >■: 

I 


(see  Reference  6,  p.  15)  we  get 
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0,  (x,  y)  in  W 


(9) 


t 


(x,  y,  z)  +  ik?  (x,  y,  z) 


]i 


0+ 

T.-0- 


Tnis  condition,  plus  Equation  {7}  applied  on  the  two  sides  of  !>,  plus  Equation 
(8),  determine  <p  as  a  solution  of  Equation  (5). 

The  conditions  satisfied  by  ?  (x,  y,  z )  are  satisfied  also  by  (x,  y,  -*). 
Hence,  Tf  is  an  odd  function  of  z-  This  implies  that  is  zero  in  the  wake. 

In  the  half  space  z  >  0,  ^  ia  a  solution  of  Equation  (5),  which  satisfies 
Equation  (8)  and  the  boundary  conditions 

?  (x,  y,  0+)  =  ~  +  ikf,  (x,  y)  in  S  (10) 

z  ox 

*x(*.  Y>  0+)  +  ikp(x,y,0+)  =  0,  (x,  y)  in  W  (11) 


■p(x,  y ,  0+ )  =  0,  (x,y)  not  in  S  +  W  (12) 

Such  a  solution  may  be  built  up  from  a  doublet  distribution  over  S  +  W  or  a 
source  distribution  over  the  half  plane  z  =  0,  x  >  0. 

4.  BASIC  SOURCE  AND  DOUBLET  SOLUTIONS  OF  THE  DIFFERENTIAL 
EQUATION  (See  References  7.  8,  and  9.  ) 

The  solution  of  Equation  (5)  which  represents  a  point  source  at  the 
origin  is 


f0.  x  S  0 


WQ  <x*  Y-*> 


2w  x 


(13) 


(See  Reference  9.  )  The  potential  of  a  point  doublet  oriented  parallel  to  the 
z-axis  is  obtained  by  differentiation,  as 


fO,  xSO 


fjfr.y.*) 


i  k  a 

2*  2 

x 


x  >  0 


(14) 


It  is  easily  verified  that  these  functions  satisfy  Equation  (5)  for  x+0.  They 
are  poorly  behaved  at  x  =  0  for  real  values  of  k. 
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To  improve  the  behavior  of  #Q  and  at  x  -  0.  aatume  that  k  has  a 
small  negative  imaginary  part.  This  causes  f0  and  ei  to  approach  zero 
exponentially  as  x  —  0+,  except  at  the  origin.  All  partial  derivatives  of  all 
orders  have  the  same  property.  Thus,  and  are  solutions  of  Equation 
(5)  everywhere  except  at  (0,  0,  0).  In  the  final  formulas  to  be  obtained, 
the  imaginary  part  of  k  can  be  put  equal  to  zero. 

Solutions  of  Equation  (5)  for  z  >  0  which  satisfy  Equation  (8)  are  given 
for  a  distribution  of  sources  as 

*gU.y,  *)  ~  f  f  n)  i>0(*-fc.  y-n.  zjdfcdq  (15) 

i>  0 


and  for  a  distribution  of  doublets  as 

*d(x,y,z)  -  JJ  p(6.n) y-,n*  *)d^dq  (16) 

i>  0 

where,  to  be  completely  general,  p(£,q)  may  be  any  function  such  that  the 
integrals  exist.  From  the  form  of  yo  and  yj,  the  region  of  integration  may 
be  restricted  to  the  plane  strip  0  <  £  <  x.  It  is  shown  in  Appendix  I  that  these 
functions  satisfy  the  following  boundary  conditions  for  z  =  0,  x  >  0: 

S  (x,  y,  0+)  =  p(x,y)  (17) 

8* 


lPd(x.  y,  0+)  =  p(x,y)  (18) 

(in  fact,  if  the  same  function  p  is  used  in  both  integrals,  =  dy^/dz). 

5.  THE  DETERMINATION  OF  f  BY  A  SOURCE  DISTRIBUTION 

One  method  of  attack  on  the  problem  of  finding  y  is  to  set  y  =  y8. 

Then,  in  terms  of  the  upwash 

w(x,y)  r  y  (x,  y,  0+)  (19) 


we  have  from  Equations  (17)  and  (15) 


f(x,  y.  z) 


fj  w (t , tj )  vn.  *)  d  £  d  t| 

i>  o 


(20) 


for  z  -  0. 
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The  values  of  w  on  S  are  known  by  Equation  (10).  Elsewhere,  w  is 
unknown,  and  it  must  be  chosen  so  that  the  boundary  conditions 
(11)  and  (12)  are  satisfied.  We  may  take  the  limit  as  z  -*  0+  in  Equation  (20) 
by  taking  the  limit  under  the  integral  sign: 

p(x,y,0+)  =  jj  w  (4,f| )pQ(x-4,  y-q,  0)d4dq  (21) 

4  >  0 

From  Equations  (11)  and  (12)  are  obtained  the  system  of  integral  equations 

jj  n)  ro(x“4»  y-n.  0)d£dq  =  0,  (x,  y)  not  in  S  +  W  (22) 

i>  o 


jj  w(4,q)  po(x-£.  y-n.  0)  d 4  d q  =  0,  (x,  y)  in  W  (23) 

4>  0 

Solution  of  Equations  (22)  and  (23),  followed  by  evaluation  of  p  according  to 
Equation  (21),  would  yield  the  values  of  p  on  S,  from  which  pressures  and 
forces  can  be  computed. 


A  box  method  based  on  a  source  distribution,  described  briefly  in 
Reference  9,  has  been  used  by  Weatherill  at  the  Boeing  Company.  Some  of 
his  preliminary  results  are  given  in  Reference  9. 

6.  THE  DETERMINATION  OF  p  BY  A  DOUBLET  DISTRIBUTION 
If  we  set  p'  -  pj,  then  by  Equations  (18),  (16),  and  (12) 

p(x,y.z)  =  f  j  ?(4.q,0+)p1U-4.y-tl,*)d4dq  (24) 

S+W 


In  terms  of 


+(x.  y) 


lim  —  p  (x,y,  z) 
z—0  *  1 


0,  x  S  0 


ik  1 

2t  2 
x 


x  >  0 


(25) 


i 
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the  normal  derivative  of  ^  at  z  =  0  is  given  by  a  singular  integral: 

=  j f  ?(£,*! .  o+)  4^<x- 4 . y-n )  d  t  d  (26) 

s+w 

The  values  of  ^  (£,  q  ,  0+)  rnust  be  determined  then  from 

ff  ?(4.  T1  °^)  (x-4,  y-T)  )  d  4  d  =  w(x,y),  (x,  y)  in  S  (27) 

S+W 

*(x,  y,0+)  =  0,  (x.y)inW  (28) 

7.  A  COMPARISON  OF  THE  METHODS 

Except  for  the  singularity  of  the  integral  in  Equation  (27),  all  points  of 
difference  are  in  favor  of  solving  the  problem  by  doublets.  There  are  these 
points: 

a.  The  region  of  integration  in  the  source  method  extends  theoreti¬ 
cally  to  i^in  t|  ;  even  practically,  the  region  must  be  extended  an 
extreme  distance.  In  the  doublet  method,  the  region  is  restricted 
to  S  +  W.  This  distinction  is  not  so  great  for  supersonic  flows. 
There,  the  region  of  influence  of  the  wing  is  swept  back  along 
Mach  lines,  and  the  set  of  points  in  this  region  that  influences  the 
wing  is  bounded  (see  Reference  10). 

b.  After  the  unknown  function  under  the  integral  sign  is  known,  the 
source  method  requires  an  extra  step  —  the  evaluation  of  ^  on  the 
wing  from  Equation  (21). 

c.  If  values  in  the  wake  must  be  considered,  the  condition  in  the 
wake  for  the  source  method.  Equation  (23),  is  more  complicated 
than  the  corresponding  condition,  Equation  (28).  lor  the  doublet 
method. 

The  doublet  method  was  used  because  of  point  a 

8.  THE  ADVANTAGE  OF  A  STRAIGHT  TRAILING  EDGE 

Suppose  the  wing  has  a  straight  trailing  edge  perpendicular  to  the 
direction  of  flow  (x  =  constant  along  the  edge):  then  the  wing  is  not  influenced 
by  the  wake.  This  is  reflected  in  the  equations  by  the  fact  that  the  inte¬ 
grands  are  eero  when  |  >  x.  Hence,  in  either  method,  for  the  determina¬ 
tion  of  ^  on  the  wing,  the  condition  in  W  need  not  be  used. 
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9.  THE  DOUBLET  BOX  METHOD 


Consider  a  flow  at  Mach  1  past  an  oscillating  wing  with  its  nose  at  the 
origin,  lying  approximately  in  the  xy- plane,  with  x  =  1  along  the  trailing 
edge.  The  value  of  the  unsteady  potential  ^  on  the  wing  may  be  found  by 
solution  of  Equation  (27),  which  may  be  written  as 

//<?  ii‘  •  0+)  4-  (x-£,  y-Tj  )  d  £  d  n  =  w(x,y),  (x,  y)  in  S  (29) 

S 

To  get  an  approximate  solution  of  this  equation,  let  the  xy- plane  be  covered 
with  a  grid  of  square  boxes  with  sides  of  length  d,  so  that  box  edges  lie 
along  the  coordinate  axes  (see  Figure  2).  Let  the  region  B  be  composed  of 
all  boxes  whose  centers  lie  in  S;  B  is  an  approximation  to  S  by  boxes.  Let 
i,  j  be  box  indexes  in  the  x-  and  y-directions.  Approximate  "Tf  by  a  constant 
value  ^  in  the  (i,  j)-th  box  B^.  Impose  the  condition  of  Equation  (7)  at  the 
center  (x^,  y j )  of  each  box  B^j  in  B,  with  the  region  of  integration  replaced 
by  B.  Then  Equation  (29)  gives  a  system  of  linear  algebraic  equations  for 
the 


// 

3., 

t'j’ 


)  d  i  d  q 


w(xityj 


(30) 


Examination  of  the  integral  in  Equation  (30)  shows  that  it  depends  on 
i,  j,  i\j’  only  via  i-i\  The  notation 

A(i-i,.)j-j*|>  =  //  ^(x.-e.yj-T,)  d  4  d  n  (31) 


is  introduced.  Formulas  for  the  evaluation  of  this  quantity  are  given  in 
Appendix  IL 


Segregating  the  terms  with  i'=i  on  the  left.  Equation  (30)  becomes 


^  A(0.  |j-j’|  )  =  wlx-.y^)-  ^  J  A(i-i\  |j-j'| )  *.tjl  (32) 

j'  i'<i  j' 

For  fixed  i  and  varying  j,  this  is  a  smaller  system  of  equations  that  may 
be  solved  for  each  consecutive  value  of  L 


11 


12 


Now  suppose  the  wing  is  symmetric  about  the  x-axis;  then  only  modes 
of  oscillation  that  are  symmetric  or  antisymmetric  in  y  need  be  treated. 
Consider  a  symmetric  mode,  will  have  the  same  value  at  corresponding 
boxes  across  the  x-axis.  This  may  be  used  to  reduce  the  range  cl  the  sumtj 
in  Equation  (32)  and  the  range  of  j.  Let  j  =  1  in  the  row  of  boxes  in  which 
0  <  y  <  d.  Then,  combining  terms  for  symmetrically  placed  boxes, 

X  [A(0.  |j-j'|)  +  A(0.j+j'-l)|  (33) 

j'il 

=  w(xi(  Yj)  -  X  X  +  J  2  1 

i'c  i  ’*2  1 

The  equations  for  j  i  0  are  implied  by  those  with  j  2  1.  Thus,  the  sine  of 
the  system  has  been  reduced  by  a  factor  of  2. 

For  antisymmetric  modes.  Equation  (33)  applies,  with  the  sums  of 
values  of  A  replaced  by  differences. 

10.  EXTENSIONS  OF  THE  METHOD 

The  computer  program  discussed  in  Section  3  has  some  restrictions 
that  are  not  inherent  in  the  box  method,  such  as  the  requirement  of  a  straight 
trailing  edge.  Some  possible  modifications  that  extend  the  applicability  of 
the  program  will  now  be  described. 

To  modify  the  program  for  modes  antisymmetric  in  y,  it  is  only 
necessary  to  change  some  of  die  signs  in  Equation  (33),  as 
indicated  in  the  discussion  above,  and  replace  even  powers  of  y 
by  odd  powers  in  the  formulas  used  for  deflection  and  potential. 

To  deal  with  a  more  general  trailing  edge,  it  is  necessary  to 
use  the  values  of  in  the  wake.  For  fixed  y,  if  x  =  xq*  at  the 
trailing  edge.  Equation  (28)  may  be  integrated  to  give 

+<x.  y,  0  +)  =  e  ‘  XT>  (xT.  y.  0  ♦) 

in  W.  In  addition  to  the  set  of  boxes  B  on  the  wing,  a  corres¬ 
ponding  set  of  boxes  Byy  on  W  must  be  considered.  After  finding 
a  value  in  a  box  of  B  along  the  trailing  edge,  the  formula 
above  may  be  used  to  find  values  in  the  boxes  directly  down¬ 
stream.  If  the  ith  row  of  boxes  includes  boxes  of  B^y,  to  the 
right  side  of  Equation  (33)  must  be  added  the  contribution  of  all 
boxes  Bi'  j/  in  Byy  with  i  <  i.  The  computer  program  must  also 
be  modified  in  several  other  respects,  to  take  into  account  the 
more  general  wing  shape. 
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A  wing  that  consists  of  several  almost  planar  sections  indifferent 
planes,  such  as  a  wing  with  folded  tips,  may  also  be  handled  by 
the  doublet  box  method.  Equation  (33)  applies,  if  ij  i®  inter¬ 
preted  as  one-half  of  the  discontinuity  in  $  between  the  upper  and 
lower  surfaces.  The  influence  coefficients  involved  are  given  by 
a  more  general  formula  (not  given  in  this  report),  allowing  for 
out-of-plane  influence  of  the  doublets.  Formulas  analogous  to 
those  ox  Appendix  II  may  be  developed,  which  are  not  much  more 
complicated.  The  main  effect  of  this  extension  on  the  computer 
program  would  be  a  greater  number  of  distinct  values  of  the 
influence  coefficients,  so  that  it  would  not  be  possible  to  store 
them  all  in  an  array  in  core  unless  the  limit  on  the  number  of 
boxes  in  each  direction  were  considerably  reduced. 

Rectangular  boxes,  not  necessarily  square,  may  be  used.  Let 
the  boxes  have  sides  of  length  dj  chordwise  and  d2  spanwise. 

If 


*  1  =  kdj ,  I  2  =  kd2/di 


the  formula  for  the  influence  coefficients.  Equation  (39)  in 
Appendix  II,  must  be  replaced  by 


A  (n,  m)  - 


iv-m |  <1/2 
lu-nl  <1/2 
u>0 


dudv  e-l/2  i  (llU  +  i2v2/ u) 
u2 


This  may  be  evaluated  by  the  methods  of  Appendix  II.  Except  for 
this  difference,  the  method  is  essentially  the  same.  The  best 
choice  of  box  shape  probably  depends  on  the  aspect  ratio  of  the 
wing. 
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3.  DESCRIPTION  Of  THE  COMPUTER  PROGRAM 


1.  COORDINATE  SYSTEMS 

An  initial  coordinate  system  (x,  y,  z)  is  assumed,  with  the  x-axis  in  the 
direction  of  the  flow.  The  undisturbed  position  of  the  wing  is  in  a  region  S 
in  the  xy-plane,  with  the  x-axis  along  the  center  line  and  the  origin  at  the 
nose  (see  Figure  1).  This  coordinate  system  is  used  in  the  data. 

In  the  program,  a  dimensionless  coordinate  system  (x,  y,  z)  is  used, 
based  on  the  root  chord  length  b: 

x  =  x/b 
y  =  y/b 
z  =  x/b 

2.  WING  GEOMETRY 

The  wing  is  symmetric,  with  trailing  edge  x  =  h.  To  complete  its 
description,  the  portion  of  the  leading  edge  on  which  y  >  0  must  be  speci¬ 
fied.  This  is  done  by  giving  the  coordinates  of  the  end  points  of  NS  line 
segments  along  the  edge  (1  -  NS  -  3},  beginning  at  a  point  at  which  y  =  0: 
(0,yo).  (xj.yj).  •  •  *  .  (xNS,  9nsV  The  edge  of  S  includes  the  polygonal  line 
through  these  points.  If  yG  >0,  it  also  includes  the  line  from  the  origin  to 
(0,yo).  If  xNS<  **  includes  the  line  from  (*ns-  yjsjs)  to  O^Yns)*  (Se® 
Figure  3.  ) 

Leading  edges  of  fairly  general  shape  may  be  approximated  by  such 
polygonal  lines. 

3.  THE  DEFLECTION  DATA 

A  mode  is  specified  by  the  vertical  deflection  function  f(x,  y)  in  terms 
of  which  the  equation  for  the  instantaneous  position  of  the  planform  is 

£  =  Re  (6  ’  eiwt  f(x,  y)J 

where  6  is  a  constant. 

In  the  program,  f  is  assumed  to  be  a  polynomial  in  x  and  y^.  The 
data  may  give  either  the  coefficients  of  this  polynomial,  or  values  of  f  at  a 


i 

J 

1 


9b* 

t 
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set  of  points  on  the  wing.  In  the  latter  case,  a  polynomial  is  fitted  to  the 
given  values  by  a  least  square  error  technique. 

4.  LEAST  SQUARE  SURFACE  FITS 

The  problem  involved  here  is  the  approximation  of  a  function  of  x  and  y 
by  an  expression  of  the  form 

7(x.y)  =  ^  anm  x“yZm  F(Xfy) 

n,  m 

when  a  set  of  values  of  the  function  is  known.  This  arises  in  the  program 
in  two  places.  The  subroutine  DREO  fits  a  representation  of  the  deflection 
of  this  type  with  F  =  1  to  the  given  deflection  values.  In  the  subroutine 
B6XP,  such  a  fit  is  made  for  the  potential,  with 


/  2  ,  .2 
v  x  -  XQ(y) 

J  1  -  y2/y2 

7  7  max 

H 

u, 

or 

-  .  i 

or 

>/x  -  Xo(y) 

1 
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(x  =  XqIy)  is  the  equation  of  the  leading  edge)  depending  on  the  wing  shape. 
This  factor  approximates  the  proper  behavior  of  p-  at  the  edges.  The  factor 

Jx*"  -  xQ(y)2  is  used  for  a  pointed  nose  (yc  =  0),  and  Jx  -  xc(y )  for  an 

unswept  nose  (yQ  >  0).  The  factor  ^1  -  y^/y^, 
has  a  side  edge  along  which  y  =  ymax. 


max 


is  included  if  the  planform 


The  factor  F(x,y)  is  real,  so  the  values  of  e  have  real  and  imaginary 
parts  that  involve  only  the  corresponding  parts  of  the  ajj^s.  Hence,  these 
real  and  imaginary  parts  may  be  handled  separately,  reducing  the  problem 
from  one  in  complex  numbers  to  (me  in  real  numbers. 

Let  arnm  =  Re[anmj,  and  let  the  real  parts  of  given  values  of  the 
function  at  data  points  be  at  (xj,  yj),  j  =  1,  ....  NP.  Then  for  the  read 
parts  we  wish  to  have 

V  a  xny  2xn  F<x.,y.)9f^.,  j  =  1 . NP 

?  nm  j  7j  J  J 

n,  m 

The  least  squares  method  minimises 

Q  = 


Z.  Z.  nm  j  'j  J  'j  J  I 

•  _  ____ 


J  n,  m 

(See  Reference  11,  Chapter  16.  ) 

For  condensed  notation,  let  r  be  a  single  index  over  the  pairs  (n,m), 
let  ohm  *  «r.  a®11  Yj2m  F(xj.  y j )  =  Ajr.  Then 


A.  - 
jr  J 


J  '  *■ 

Let  the  range  of  r  be  from  1  to  NC  -  NP. 
To  minimise  Q,  we  set 


8Q 

8a 


=  0,  r  =  1 . NC 


i 
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This  leads  to  the  system  of  equations 


Put 


Then  Equation  (34)  reduces  to 


2 


NC 


(34) 


(35) 


(36) 


The  matrices  (Brr.)  and  (C,)  must  be  set  up  to  solve  Equation  (36). 

It  is  not  necessary,  however,  to  set  up  the  matrix  (Ajr).  Only  one  row  of 
(Ajr)  is  needed  at  a  time.  This  is  fortunate,  because  the  program  allows 
(Ajr)  to  become  as  large  as  2500  x  20.  For  each  value  of  j,  the  jth  row  of 
(Ajr)  is  computed,  and  from  this  the  jth  terms  in  the  sums  in  Equation  (35) 
are  formed  and  added  in. 


In  the  complex  case,  there  is  a  corresponding  system  of  equations  for 
the  imaginary  parts: 


2 


r' 


B 


rr 


C" 

r 


The  two  systems  of  equations  are  solved  together  by  the  subroutine  XS1MEQ, 
which  allows  for  more  than  one  set  of  values  on  the  right. 

5.  GENERALIZED  FORCES 


The  generalized  force  coefficient  is  defined  (Reference  6)  by 


L. . 

U 


1 

1/2 pU  2S 

S3 


//ap  ^x,  y)f.(x,y)dxdy 
S 
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where  APi  is  the  lifting  pressure  difference  in  the  ith  mode,  and  fj  is  the 
deflection  function  in  the  jth  mode.  In  terms  of  the  potential  ?(x,y)  on  the 
upper  surface, 


Ap.  =  2pU*  (*x  4  i  k  r) 

Ly =  s  //'% Hkf,,id‘dy 
s 


After  integration  by  parts, 


L. . 

‘J 


_4 

S 


/  *fjd'r+  //*(ikfj-'»^)d,‘<iy 
x=l  S 


(37) 


In  Equation  (37)  insert  the  series 


—  -  V  n  2m  . 

*  -  2  *»-  y  F(x'y> 


n,  m 


*i-  I 


n*  2m' 


1  .  .  *  y 

n'm'  ’ 


n' ,  m' 


The  result  is 


=  |  y  d ,  ,  y 

ij  S  ZL  n'm'  Z, 


nm 


n',  m' 


n,  m 


1  f  2m+2m' 

t  y 


2  ; 

X  =  1 


F(l,  y)  dy 


♦  ik 


i  // 


n+n1  2m+2m'  ,  .  . 

x  y  F(x,  yjdxdy 


.  1  f  f  n+n'- 1  2m+2m' 

-  n'  •  —  x  y 


F(x,  y)  d  x  d  y 


The  integrals  in  this  expression  depend  only  on  the  wing  shape.  They 
are  computed  by  the  subroutine  F0RCI  before  the  work  on  the  individual 
modes  begins.  During  the  work  on  the  ith  mode,  the  sum  over  n  and  m  Is 
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performed  in  the  last  part  of  the  subroutine  B^XP,  for  each  set  of  values  of 
n*  and  m'.  The  sum  over  n'  and  m'  and  multiplication  by  8/S  is  performed 
in  the  last  part  of  the  main  program. 

6.  THE  USE  OF  GAUSSIAN  QUADRATURE  IN  THE  EVALUATION  OF 
GENERAUZED  FORCES 


Gaussian  quadrature  is  an  approximation  of  the  form 


N 


/ 


f(u)  du  3T 


V'V 


j=l 


exact  for  polynomials  of  degree  -  2N  -  1.  (See  Reference  11,  Chapter  7.  ) 
This  formula  is  used  with  (a,  b)  =  (0,  1 ),  N  =  6.  The  values  of  the  hj's  and 
uj's  for  this  case  were  obtained  from  values  listed  in  Reference  11.  They 
are  given  as  H(l) . H(6),  U(l), ....  U  (6)  in  the  subroutine  SECT. 

Subroutine  F^RCI  finds  the  values  of 


AXY(I.J)  =  j  fj x1'  1  y2J'  2  F(x.y)  d  x  d  y 


AY  (J )  =  J  /  y2J  2  F(1 ,  y)  d  y 
x=l 

for  I,  J  =  1,.  .  .  ,  9.  To  do  this,  the  contributions  to  the  integrals  from  each 
section  of  wing  behind  a  straight  piece  of  leading  edge  are  calculated 
separately  in  SECT. 


The  form  of  F(x,  v)  is 


F(x,  y) 


r 

/x  -  xQ(y) 

yr~7v 

'  '  max 

or 

.  .  , 

or 

/  2  .2 
Vx  -  x  (y) 

° 

1 

*  ^ 

depending  on  the  wing  shape.  We  have  integrals  that  behave  like  square 
roots  at  the  leading  edge.  The  integrals  over  one  wing  section  are  of  the 
form 
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y+ 


BXY 


/  ■*»  / 

y-  »0(y) 


I-  1  2J-2 

ax  X  y  F(x,y) 


And 


y+ 


2J  -  2 


BY  *  J  dyy““  “  F(l,  y) 

y_ 


In  BXY,  the  chordwise  integral  is  evaluated  first  at  each  value  of  y 
at  which  it  will  be  needed.  The  new  variable 


=  A  -  *o(y) y'  yi  -  *0<r) 


(38) 


is  introduced.  Then 

j  dx  *  y^J  2  F(».  y)  *  y  *1-2^1  -  *o(y)l  x1  1  y2j  2F<x.y) 
xo(y»  0 

The  integrand,  as  a  function  of  u,  is  well-behaved  at  the  leading  edge.  It 
is  approximated  by 


6 

e(y)  =  J  hi*  2 [* "  *0w]*i1' 1  y2J  F(*i»y) 

i«l 


where  x^  is  computed  from  the  value  of  u^  according  to  Equation  (38), 

In  the  y- integration  in  BXY  and  BY,  the  integrand  approaches  sero  as 
as  y  —  Yxrax  lik*  -  y/ymax  or  (1  -  y/Tmax^2*  Accordingly,  the  change 
of  variable 


y  = 


y+-«y+-y>.  T,<yMX 
y.  -  <y+  -  yV-  n  *  W 


is  used,  which  makes  the  interval  of  integration  0  <  v  <  1  and  removes  the 
square  root  behavior  in  the  last  section  of  the  wing.  This  leads  to  the 
formulas 
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BXY 


'  (yt  -  y.)  2.  hj  "  *y j  >  ‘ 


'•  y+<  ymax 

2Y  y+  =  y 


max 


BY  =  (y+  -  y.>  ^  hj  v/J"2  F(l.y.)- 
j=l 


*•  y+  <  ymax 

2u..  y+  =  y 


max 


7.  LEADING  EDGE  CORRECTION 


The  value  of  potential  found  for  each  box  from  Equation  (33)  is  taken 
to  be  the  value  of  ^  at  the  box  center.  Thus,  the  values  obtained  are  in 
error  only  by  virtue  of  the  error  introduced  in  the  values  of  upwash  when 
the  actual  distribution  of  potential  in  a  box  is  replaced  by  this  constant 
value.  This  error  is  especially  important  in  the  first  row  of  boxes,  for  a 
wing  with  an  unswept  leading  edge.  The  major  effect  is  on  the  upwash 
values  in  that  row. 


To  estimate  this  error,  consider  the  two-dimensional  case,  in  which 
?  is  independent  of  y.  In  Equation  (26),  the  expression  for  upwash  due  to  a 
doublet  distribution,  integrate  by  parts  over  £,  then  integrate  over  q.  The 
result  is 


w(x.y)  =  ““//  0+)  c 


0<4<  x 


(x  -  !)* 


--ikfx  -  t  +  1* 

1  ff  dtdn  /__  .  1  ..  2  \  6  *’i  I 

=  7  JJ  -“2K+2lk')C 


C<4<* 


<y  -  ) 


.ylHK  T  -■ 

V  n  J  \/*  ~  4 


-Tik(x  -  4) 


(q  +ilk*) 


For  x  *  yd,  if  kd  is  small, 
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*£ (i.n.  o+) 


The  correct  leading  edge  behavior  is  possessed  by  the  expression 
y  -  C  .  We  have 


If  ?  is  constant  on  0  <  £  <  d,  and  has  the  value  C  \f\  !Z  d,  then  in  the 
above  integral,  can  be  expressed  in  terms  of  a  delta  function: 


Accordingly, 


Note  that  the  latter  value  is  smaller  than  the  value  of  "w  evaluated  for 

~  Cjj by  the  factor  2/ir.  This  implies  that  the  values  of  potential  found 
for  the  first  row  of  boxes  will  be  more  accurate  if  the  upw ashes  in  that  row 
are  multiplied  by  2/*r. 


8.  THE  FORM  OF  OUTPUT 


The  viewpoint  is  taken  that  calculation  o 4  the  generalised  forces  is  the 
basic  purpose  of  the  program.  They  are  always  printed  out.  There  are 
other  outputs  that  will  be  printed  if  the  appropriate  data  signal  is  given. 

Each  of  the  following  is  printed  if  the  data  item  specified  in  parentheses  is 
non-sero: 


a.  The  coefficients  of  the  deflection  polynomial,  if  it  has  been  com¬ 
puted  as  a  fit  to  given  values  of  deflection,  DA(87) 

b.  The  upwash  array,  DA(88) 

c.  The  potential  array,  DA(89) 

d.  The  coefficients  of  the  potential  series,  DA(90) 
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e.  Values  of  pressure  and  potential  at  the  box  centers,  computed 
from  the  series,  DA{91). 

9.  THE  DATA  SUBROUTINE  DATRD 

This  subroutine  reads  all  data  items  into  the  array  DA.  Punched  cards 
used  for  data  are  considered  to  contain  six  fields  of  length  12  as  indicated  in 
the  sample  data  sheets.  The  first  field  contains  information  for  DATRD. 
Ending  in  column  12  is  an  integer  giving  the  location  in  the  data  array  for 
the  entry  in  the  second  field.  The  following  fields  go  into  consecutive  loca¬ 
tions,  if  the  data  are  numeric.  Floating  point  numbers  should  be  written 
with  decimal  points,  and  fixed  point  numbers  adjusted  to  the  right  end  of  the 
field. 


The  word  ALPHA  in  columns  2  through  6  indicates  that  the  data  on  the 
card  are  alphanumeric.  These  are  stored  in  DA  in  a  different  way,  taking 
up  ten  locations  per  card.  The  data  may  be  printed  later,  just  as  they 
appeared  on  the  card. 

On  a  numeric  card,  if  a  field  is  blank  the  corresponding  location  in  DA 
is  unchanged.  This  is  not  true  for  an  ALPHA  card. 

A  minus  sign  in  column  1  indicates  the  last  card  to  be  read  at  the  time. 
DATRD  reads  cards  until  this  minus  sign  is  encountered,  then  returns  to  the 
main  prog  ran:. 

10.  A  NOTE  ON  THE  USE  OF  TAPES 

In  writing  of  tins  program,  the  following  tape  numbers  have  been  used: 
output  tape,  number  6;  input  tape,  number  5;  and  tape  simulated  by  an  internal 
file,  number  99. 

The  tape  numbers  5  and  99  appear  in  the  subroutine  DATRD.  Elsewhere 
only  the  output  tape  number  is  used.  It  occurs  in  the  main  program,  and  in 
the  subroutines  SHAPE,  DRED,  B0XP  and  B@XP($. 

11.  USE  OF  THE  PROGRAM  FOR  FIXED  WING  AND  MODES  AT 

VARIOUS  FREQUENCIES 

If  a  non-zero  quantity  is  entered  in  the  appropriate  location  in  the  data 
array,  (DA26),  it  indicates  that  a  wing  shape  and  set  of  modes  to  be  used 
are  the  same  already  used  for  another  frequency.  Then  quantities  that 
depend  only  on  wing  shape  and  deflection  data  will  not  be  computed,  but  will 
be  taken  from  the  permanent  arrays  in  which  they  were  stored  in  the  pre¬ 
vious  case.  The  number  of  boxes  along  the  root  chord,  DA(27),  may  not  be 
changed  when  this  is  done. 
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When  this  option  ie  exercised,  ail  work  for  the  present  frequency  will 
be  carried  out  after  reading  one  set  of  data,  which  need  only  include  the 
frequency  and  the  indicator,  DA  (26X  Titles  for  the  individual  modes  are 
not  printed. 

12.  DESCRIPTION  OF  THE  DATA  ARRAY 

All  data  are  entered  into  the  array  DA,  dimensioned  for  700,  as 
described  in  Paragraph  9.  Th'n  layout  of  the  array  is  as  follows: 


1  -  10 

Title 

13  -  2c 

Mode  Title 

23  : 

Frequency  (cycles  per  unit  of  time),  v 

24  : 

Root  chord  length,  b 

25  : 

Speed  of  sound,  a 

26  : 

Indicator  for  new  frequency  (See  Paragraph  11) 

27  : 

Number  of  boxes  along  root  chord 

28  : 

Number  of  modes 

29  : 

Number  of  sections  of  leading  edge  to  be  given 
(See  Paragraph  2) 

30  -  36 

Coordinates  of  points  on  the  leading  edge 
(See  Paragraph  2) 

39  : 

Indicator  to  suppress  calcination  of  potential  for 
a  mode 

46  -  70 

Coefficients  of  the  deflection  polynomial 
(See  Paragraph  3) 

87  -  91 

Output  indicators  (See  Paragraph  8) 

98  : 

Number  of  points  at  which  deflections  are  given 
(See  Paragraph  3) 

99  : 

Number  of  x  values 

100  : 

Number  of  y  values 

101  -  700 

Deflection  data  for  a  maximum  of  150  points 

Note:  23,  24.  25  end  30-36  must  be  entered  in  consistent  unite 
of  length  end  time. 

13.  OUTLINE  OF  THE  PROGRAM 

For  the  purpose  of  description,  the  main  program  has  been  divided 
into  20  parts,  as  indicated  in  Figure  4,  which  shows  the  flow  of  the  program 
end  the  subroutines  called. 


14.  SIZE  LIMITATIONS  OF  THE  PROGRAM 

The  program's  size  limitations  are  as  follows: 

a.  Box  size  —  the  half  wing  must  be  enclosed  in  a  rectangle  that 
contains  no  more  than  50  boxes  in  each  direction.  (The  use  of  a 
large  number  of  boxes  is  not  recommended,  because  the  time 
required  is  roughly  proportional  to  the  cube  of  the  number  of 
boxes  along  the  root  chord.  The  possibility  of  50  boxes  in  each 
direction  is  intended  to  allow  a  large  range  in  aspect  ratio.  ) 

b.  Number  of  modes  —  ten  at  most. 

c.  Number  of  points  at  which  the  deflections  are  given  for  one 
mode  —  150  at  most. 

d.  Terms  in  the  deflection  polynomial  —  this  is  5Zd  xny^mt  where 

„  ,  .  nm  7 

0  s  n  <  4,  0  s  m  s  4. 
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I3> 


DATRD 


1  SHAPE] 


^  Fjrwci 

fp0T2l 

© 

© 

^  DATRD 

]  © 

© 

© 

•1  dred| 

© 

© 

© 

►}wval| 

© 

© 

b^xp|- 

♦♦j*0XP0|(T8) 

(7)  The  doto  orroy  is  initiolized  by  making  the  titles  blank  and  weij^h  1  .0. 

^2^  Subroutine  DATRD  reods  the  first  block  of  data 

^3^  The  following  quantifier  ore  computed:  reduced  frequency 
k  -  2  w vb/a,  D  =  dimension lers  box  size,  DH  3  0.50. 

(4  J  The  title  is  printed,  with  the  volues  of  L  (the  number  cf  boxes  along  the  chord), 
k,  and  v. 

^5^  If  DA(26)  i  0  (See  Paragraph  1]  of  this  section)  control  transfers  to  Part 

Subroutine  SHAPE  computes  the  number  of  boxes  in  eoch  row  and  the  wing  area. 

Subroutine  F0RCI  computes  the  integrals  to  be  used  in  evaluation  of  generalized 
forces  (See  Porogropht  5,  6  of  this  section). 


DATRD  reads  the  next  block  of  data. 


given.  The  <  oefficients  ore  transformed  into  coefficients  of  o  polynomlnol  In  the 
dimensionless  variables  x,  y  and  stored  In  the  permanent  orroy  DF. 


©  G(M)  -  DA(37) 


current  mode.  In  this  cose,  this  mode  will  be  used  only  a  o^eflection  In  the 
computation  of  geiaralized  forces.  Control  goes  to  Part  fl9) 

Subrout  ine  WVAL  sets  up  the  array  of  downwash  values  W. 


of  this  section. 


the  orroy  of  potential  values.  A  least-squares  fit  is  made  to  these  volues  ns 
described  in  Parograph  8  of  this  section.  The  print  options  (b),  (c),  (d),  (e)  listed 
in  Paragraph  8  of  this  section  are  handled  by  this  subroutine.  The  final  section  of 
the  subroutine  stores  in  the  orroy  PS  the  Sims  over  tfse  potential  coefficients  a^ 
drown  in  Porogroph  5  of  this  section. 


©» 


there  are  more  modes,  control  returns 


to  Part  (7) 


The  generalized  forces  are  computed  by  summing  aver  the  deflection  coefficients 
at  detcribad  in  Paragraph  5  of  tbit  wction.  Control  poo  to  Port  f2) 


Figure  4.  Program  Flow  Diagram 
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4.  RESULTS 


1.  THE  ASPECT  RATIO  1.  5  DELTA  WING 

The  computer  program  was  run  for  the  plunging  and  pitching  modes 
(pitch  axis  at  x  =  0)  at  the  reduced  frequencies  k  =  0.  2,  0.  5,  0.  8,  1. 0.  Forty 
boxes  along  the  root  chord  were  used,  which  leads  to  about  300  boxes  on  the 
half- wing. 

Theoretical  values  for  comparison  were  calculated  from  Davies'  formu¬ 
las  (Reference  12).  These  are  analytic  expressions  of  the  solution  of  Equa¬ 
tion  (5)  for  the  potential  and  generalized  forces  for  the  delta  wing  in  rigid 
modes  of  oscillation,  expressed  as  series  in  k.  Figures  5  through  7  show 
the  values  of  generalized  forces  Ln  (lift  due  to  plunge),  L21  (lift  due  to 
pitch),  and  L22  (moment  due  to  Ditch).  Note  that  the  vertical  scales  have 
been  expanded  in  the  portion  of  interest,  especially  for  L}  |.  Most  of  the 
values  agree  to  within  2  or  3  percent. 

The  differences  indicate  the  errors  introduced  by  the  box  method  in 
the  solution  of  Equation  (5),  as  distinguished  from  the  errors  inherent  in 
this  equation. 

Figure  8  gives  the  chordwise  distribution  of  values  of  4  for  the  plunging 
mode  at  k  =  0.  5,  for  y  =  yroax^  =  125. 

2.  THE  ASPECT  RATIO  2.  0  RECTANGULAR  WING 

The  plunging  and  pitching  modes  were  again  used  at  k  =  0.  3,  0.  6,  0.  9. 
Twenty-five  boxes  were  allowed  along  the  chord,  giving  625  boxes  on  the 
half -wing.  The  values  of  L21  and  L22  are  shown  in  Figures  9  and  10,  with 
values  from  Land  ah  1  (Reference  6,  page  84)  for  comparison.  Landahl's 
values  were  obtained  by  a  method  of  solution  of  Equation  (5)  which  applies 
onl)  to  a  rectangular  wing  in  modes  of  oscillation  with  a  deflection  independent 
of  y. 

3.  THE  ASPECT  RATIO  3.  0  RECTANGULAR  WING 

Finally,  for  the  aspect  ratio  3.  0  rectangular  wing,  a  comparison  is 
made  with  experimental  pressure  values.  These  values  were  given  in 
Reference  13  for  a  5-percent  thickness  wing  oscillating  in  an  elastic  bending 
mode.  At  Mach  1,  the  reduced  frequency  was  0.24.  The  chordwise  pressure 
distribution  at  y  -  ymax^  shown  in  Figure  11. 
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Figure  6.  Lift  Due  to  Pitch  for  an  Aspect  Ratio  1.  5  Delta  Wing 
(Compared  With  Reference  12) 
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Figure  7.  Moment  Due  to  Pitch  for  an  Aspect  Ratio  1.  5  Delta  Wing 

(Compared  With  Reference  12) 
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Figure  I?.  Chordwise  Pressure  Distribution  on  an  Aspect  Ratio  3.0 
Rectangular  Wing  in  an  Elastic  Mode;  y  -  ymax/2  =  0.  75 
(Compared  With  Reference  13) 


The  variation  of  the  experimental  values  from  the  computed  values  is 
of  the  type  that  thickness  effects  should  be  expected  to  cause:  the  measured 
pressure  is  (1)  smaller  near  the  leading  edge,  (2)  larger  before  the  point  of 
maximum  thickness  (x  -  0.  5),  and  (3)  smaller  beyond  this  point.  The  experi¬ 
mentally  determined  values  of  local  Ms  h  number  along  this  chord  range 
from  0.  84  to  1.  35,  which  indicates  tha'  he  thickness  has  a  considerable 
effect  on  the  flow. 

The  theoretical  curve  given  in  k  ference  13  was  obtained  from  the  sub¬ 
sonic  kernel  function  method,  applied  at  M  =  0.  99.  This  curve  is  included  to 
show  how  another  theoretical  method  compares  with  the  experimental  values. 

4.  COMPUTER  RUNNING  TIME 

The  results  described  in  this  section  required  about  20  minutes  total 
computer  time  on  the  IBM  7094.  With  nonessential  output  omitted,  this  time 
could  have  been  reduced.  All  optional  output  was  given,  resulting  in  about 
40,  000  lines  of  output. 


5.  CONCLUSIONS 


A  procedure  has  been  developed  for  predicting  unsteady  aerodynamic 
forces  and  pressures  on  an  oscillating  wing  by  the  use  of  the  transonic  box 
method.  The  results  obtained  by  this  method  agree  quite  well  with  theoreti¬ 
cal  values  from  other  methods  that  are  applicable  only  to  special  planiorms 
The  box  method  has  the  advantage  of  applicability  to  a  general  planform. 

The  only  other  method  of  this  generality  at  Mach  1  is  the  sonic  limit  of  the 
subsonic  kernel  function  method  (see  Reference  16)  that  has  not  been  very 
successful. 

The  comparison  with  experimental  values  in  Figure  11  indicates  that 
the  most  serious  limitation  of  the  method  is  that  thickness  is  neglected. 
Thickness  may  be  incorporated  into  a  box  program  by  using  modified  forms 
for  sources  and  doublets,  depending  on  the  local  Mach  number  (see 
Reference  14).  This  was  not  accomplished  under  the  present  program. 

Other  possible  extensions  of  the  transonic  box  method,  that  would  not  require 
much  change  in  the  existing  computer  program,  are  described  in  Paragraph  10 
of  Section  2. 
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APPENDIX  I.  PROPERTIES  OF  SOURCE  AND  DOUBLET  DISTRIBUTIONS 


1.  BOUNDARY  BEHAVIOR  OF  A  DOUBLET  DISTRIBUTION 
We  wish  to  evaluate 


*d(x,y,  0+)  =  lim 
2—0+ 


The  grand  is  zero  for  £  >  x  . 


*  (x,y,  0+)  =  lim 

E”*0+ 


JJ  d£  dq  p(£.q)  *  j(x-£.  y-q.  z). 

£>0 

If  we  define  p(£,T>)  =  0  for  £  <  0,  then 

JJ  d£dq  p(4»ti)?1(x-£.  y-n.  *)  • 

£<x 


Put 


£  =  x  -  z  u 


t)  -  y  -  zv 

'hen,  using  Equation  (14),  we  have 

2 


*  (x,y,  0+)  =  lim  ff  dudv  p(x-z6u,  y-zv) ■  z3  *  (z2u,  *v,  z) 
a  z— 0+  1 


z—  0+ 


u>0 


lim  jf  dudv  p(x-z^u,  y-zv)  •  —  e 


1  2  l+v^\ 

i  xH* a*—) 


2—0+  u>0 


2«  2 

u 


Let  (x,  y)  be  a  point  at  which  p  is  continuous.  Then  the  value  of  p  in  the 
integrand  approaches  p(x,y)  as  z—  0.  Taking  the  limit  under  the  integral 
sign, 


7d(x,y,  0+)  =  p(x,  y)*^-  JJ  dudv 

u>0 


-iiki-ti 

2  u 
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Let  v 


s 


/u  .  Then 


*d<*.  y,  0+) 


0 


du 


3/2 

u 


-jik/u 


These  integrals  may  be  reduced  to  a  standard  form  by  rotating  the  paths 
integration  in  the  complex  plane  into  positions  in  which  the  exponents  are 
negative,  then  making  the  substitutions 


i  k 


The  result  is 

1 

CD  _ 

?d(x.  y.  0+)  =  p(x,  y)  •  —  f  e  P  P  2  dp 

0 


®  2 
/  e  q  dq 


=  p(x,  y) 

since  both  integrals  have  the  value  y/ir-  (see  Reference  12,  formulas  507, 
512).  Consequently,  Equation  (18)  is  valid  at  any  point  of  continuity  of 
p(x,y\ 


2.  BOUNDARY  BEHAVIOR  OF  A  SOURCE  DISTRIBUTION 


To  evaluate  *as(x,  y»  0+),  note  that 

17  =  //  dtd’'  y‘i>  •> 

t>e 

=  f I  p(£»q)  y-q*  *) 

l>c 


Hence,  by  the  result  of  the  preceding  section,  if  (x,  y)  is  a  point  at  which 
p(x,  y)  is  continuous. 


0+)  =  *d(x,y.  0+)  =  p(x,  y) 


which  verifies  Equation  (17). 


tteHNBSMltlt  mmum*  *#»■**» 


APPENDIX  II.  EXPRESSIONS  FOR  THE  INFLUENCE  COEFFICIENTS 


Equation  (31)  may  be  expressed  more  conveniently  in  terms  of 

u  =  -  £)/d 

v  =  (y.  -  q)/d 

m  =  jj  -  j*| 

n  =  i  -  i' 
i  =  kd 

(d  =  box  side  length).  By  Equation  (14)  we  have 


A(n,  m)  -  // 

|u-n|<I 
u>0 

It  is  assumed  that  /  is  small.  If  t  <  0.  1,  the  following  approximation 
gives  an  error  of  less  than  0.  1  percent  in  the  value  of  A: 

-i-ilu  "  ifn  ii (u  -  n) 
e  *  e  e 

1  -~lf (u  -  n)-i/2(u  -  n)2] 


du  dv 


4- (-4) 


(39) 


a  e 


This  reduces  Equation  (39)  to 


A(n,  m) 


-jUnj/  i  i  ZZ\ffiu 


dv 
Z~ * 


(4“4A)//^ 


-  li  v2/u 
dv  2 
—  e 


i*7/ 


du  dv  e 


-iuvJ/u 


where  the  limits  of  Integration  are  the  same  as  in  Equation  (39). 

The  following  formula  expresses  these  double  integrals  in  terms  of 
single  Integrals: 


u  v  4u 

/  d“  /  dV~P  * 

U  V 

1  1 


1  4.  2/  2 

u,  ,  --  ii  v  /u 

,  2  vdu  2 


1  f  2  vdu 
-  2p  /  uP  8 


'v  *  v. 


+  3  -  2p  / 


v  -jlfv2/u  2 

■  2  dv  2 


u  *  u, 


Equation  (40)  becomes 

,  f  A  (n  ♦  m)  -  A  (n  -  4  ,  m),  n>0 

ik  "I 

A(n,  m)  *  -r—  e  •  .  (41 

Zlr  A  (4* .  m),  n  ■  0 

!  O  £ 


48 


where 


A  (u,m) 
n 


u 


j 


du  e 


~  iiv^/u 


8 


I 


+ 


12 


<2u 


-yifv2/u 

dv  e 


(42) 


-  B  (u,  m)  +  C  (u,m) 
n  n 

Bn  and  CQ  denote  the  contributions  of  the  terms  containing  the  u  -  integrals 
and  v  -  integrals,  respectively. 

Bn(u,m)  may  be  expressed  in  terms  of  the  sine  and  cosine  integrals 

S(x)  =  /  dt 

l 

C<x)  =  /  ^^dt 
1 

(C(x)  and  S(x)  are  evaluated  by  the  subroutine  CIN .  )  The  resulting  formula 
for  Bn  is 


49 


B  (u ,  m } 

Q 


[v  (  t11  +i|2n)4 «  u3v3|  [c(^r) 


-y  ifv*/u 


(43) 


i  S 


(#)] 


m  + 


v  =  m  - 


To  evaluate  Cn(u,m),  put  v  =  8  +  m.  Expanding  part  of  the  exponen¬ 
tial  gives  the  approximation 


/ 


m+j  iivZ/u 

Z  dv  e  Z 


1 

mT 


■/. 


-4*  -4  if  (rr4  +  2ms  +  s^)/u 

c  ,  2 

ds  e 


s  e 


1  Zj  l 
-  — llm  /u  -r- 

2  f  2 


f"  dg« 
2 


-  if  ms/u 


(- 


and  performing  the  integration  gives 
1 


/ 


m- 


_  -yifv^/u  -y  if  m^/u 

W  -  C,  €rn 

dv  e  =  e 


1 

mT 


sin  ffm/Zu)  / 

1  .i  Al\ 

lm/  2u  \ 

8  uj 

1U 


f  m 


8 in  (fm/2u)  / ,  . 

- — 1 — 7— - *-  -  cos  (f m/2u) 

lm/  2u 


,  m  *  0 


For  small  valuee  of  im/2u,  the  trigonometric  functions  of  this  argument 
are  expanded  in  power  series.  To  sufficient  accuracy. 


m 


/ 


m-- 


+4  -4  if  v^/u  '4  if  m ^ / u 


2  .  2 
dv  e 


=  e 


1  (tra  \  J_  if 
6  ^  2u  y  24  u 


,  im/2u  <  0.  2 


50 


!t  may  be  verified  that  thit  it  valid  for  m  *  0. 


Combining  these  expressions. 


C  (u,m) 
n 


=  e 


-j  ifm^/u 


f  ('  4U“  i,v)+  2(t"  +i,2n)  •  tt'2u 

(‘  •  i  -  «-*.-/*-»]  . 


im/  2u  >  0.  2 


(44) 


,  1  ii  1  /im\  ,, 

1  -  —  —  -  7-  (— — I  ,  im/2u  <  0.  2 
24  u  6  \  2u  / 


The  subroutine 
Equations  (41),  (42), 


P<JT2  evaluates  tne  influence  coefficients  according  to 
(43),  and  (44). 
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COMPLEX  SIMULTANEOUS  EQUATION  SUBROUTINE 
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APPENDIX  JV.  SAMPLE  DATA  SHEETS 


The  following  pages  are  sample  data  heets  for  a  computer  run  on 
three  modes  at  three  frequencies.  The  1  itential  will  not  be  computed  for  the 
first  mode.  The  generalized  forces  fout  will  be  L<21*  L>22*  L»23»  L)l,  L»32» 

l33- 


Of  the  first  fourteen  cards,  the  cards  numbered  6,  9,  10,  11,  12,  13, 
14  do  nothing  and  are  included  only  to  indicate  how  all  data  is  entered. 

Cards  1  through  14  are  complete  in  this  respect,  and  all  later  cards  are  of 
the  same  type  as  one  of  the  first  fourteen.  The  data  used  in  the  least 
squares  surface  fit  for  the  deflection  is  entered  in  locations  101  through  J 00. 

Card  number  22  represents  56  cards  for  the  intermediate  data  points 
which  would  have  to  be  included  in  an  actual  run. 
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NUMBER  IDENTIFICATION  DESCRIPTION  DO  NOT  KEY  PUNCH 


NUMBER  IDENTIFICATION  DESCRIPTION  00  NOT  KEY  PUNCH 
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